library(knitr)
knitr::opts_chunk$set(cache = TRUE, cache.lazy = FALSE, warning=FALSE,
fig.width=6, fig.height=6, autodep=TRUE)
buildFault = function(l=400, setoff=0, rundir="TPV15", run=1) {
# through fault
x0 = 0
y0 = 0
z0 = 0
L = 28e3
W = 15e3
w = l
nl = round(L/l)
nw = round(W/w)
strike = 0
dip = 90
rake = 0
ddot = 0
sectNum = 1
system2("/Users/dinger/RSQSim.github.master/BuildFaultModels/singlePlanarFault",
args=c(x0, y0, z0, L, W, nl, nw, strike, dip, rake, ddot, sectNum),
stdout=file.path(rundir, "fault1.flt"))
# branch fault
strike = 30
x0 = setoff * sin(strike*pi/180)
y0 = 16e3 + setoff * cos(strike*pi/180)
z0 = 0
L = 12e3
nl = round(L/l)
nw = round(W/w)
strike = 30
sectNum = 2
system2("/Users/dinger/RSQSim.github.master/BuildFaultModels/singlePlanarFault",
args=c(x0, y0, z0, L, W, nl, nw, strike, dip, rake, ddot, sectNum),
stdout=file.path(rundir, "fault2.flt"))
system2("cat", args=c(file.path(rundir, "fault1.flt"),
file.path(rundir, "fault2.flt")),
stdout=file.path(rundir, paste0("faults.flt")))
unlink(c(file.path(rundir, "fault1.flt"),
file.path(rundir, "fault2.flt")))
flt = readFault(file.path(rundir, "faults.flt"))
return(flt)
}
Note that the 3 km x 3 km nucleation region centered on \(y\) = 8 km, \(z\) = -7.5 km may or may not fit exactly in a mesh with a given side length.
initStresses = function(flt, sigma=rep(120, length(unique(flt$iseg))),
tau=rep(70, length(unique(flt$iseg))),
ynuc=8000, znuc=-7500, wnuc=1500, segnuc=1, taunuc=81.6) {
sigma.out = sigma[flt$iseg]
write(sigma.out, file=file.path(rundir, "initSigma.txt"), ncol=1)
tau.out = tau[flt$iseg]
nuc = which(flt$iseg == segnuc & abs(flt$z - znuc) <= wnuc & abs(flt$y - ynuc) <= wnuc)
tau.out[nuc] = taunuc
write(tau.out, file=file.path(rundir, "initTau.txt"), ncol=1)
}
These are simple and the same for all the runs we’ll be doing here.
TPV15 (and TPV15) use \(\rho = 2670 \mathrm{kg}/\mathrm{m}^3\), \(V_s\) = 3464 m/s, and \(V_p\) = 6000 m/s. RSQSim needs the Lame parameters \(\lambda\) and \(\mu\). In terms of the density and wave speeds these are
\[\begin{align} \mu &= \rho V_s^2 \\ \lambda &= \rho \left( V_p^2 - 2V_s^2 \right) \end{align}\]rho = 2670
Vs = 3464
Vp = 6000
lameMu = rho * Vs*Vs / 1e6 # in MPa
lameLambda = rho * (Vp*Vp - 2*Vs*Vs) / 1e6
We will use our usual fixed values of \(V^\mathrm{EQ}\) (though will definitely want to try variable \(V^\mathrm{EQ}\) eventually) of 1 m/s (update: given the specified \(V_s\), \(\Delta\tau\), and shear modulus, the shear impedance relation implies a slip speed of 1.57 m/s, but I have adjusted it to 1.2 m/s to better match the rupture propagation speed seen in the dynamic models), \(V^\ast\) of \(10^{-6}\) m/s, \(D_c\) of \(10^{-5}\) m, \(a\) of 0.01, and \(b\) of 0.015. Then we will choose \(\mu_0\), and \(\theta_0\) from the two equations in http://rsqsim.ucr.edu/wiki/index.php/Stress_drops_in_artificially_nucleated_events#Further_attempt_at_more_accurate_estimate_of:
\[\begin{align} \mu_0 &= \frac{\tau_f}{\sigma} + \left(b-a\right)\ln\left(\frac{V^\mathrm{EQ}}{V^\ast}\right) \\ \theta_0 &= \left(\frac{D_c}{V^\ast}\right)^\frac{b-a}{b-a+a^\prime} \left(\frac{l b \sigma}{G V^\mathrm{EQ}}\right)^\frac{a^\prime}{b-a+a^\prime} \exp\left(\frac{\tau_p/\sigma - \mu_0}{b -a +a^\prime}\right) \end{align}\]And then force the nucleation region to fail immediately at its initial shear stress level by lowering \(\theta\) there by the appropriate amount. RSQSim can do that automatically, but I’ll just calculate the numbers here and write out an initial \(\theta\) file instead.
The approx for \(\tau_f\) is known to be pretty rough, so after an initial run showed that the final \(\tau\) was about 2.6 MPa too high, I adjusted the \(\tau_f\) in the formula for \(\mu_0\) above by this amount (in code block below).
For TPV15 (and TPV15) we have that \(\tau_p = \mu_s \sigma = 0.677 \cdot 120 \mbox{ MPa} = 81.24 \mbox{ MPa}\) and \(\tau_f = \mu_d \sigma = 0.525 \cdot 120 \mbox{ MPa} = 63 \mbox{ MPa}\).
fricTheta0 = function(flt, rundir, mus=0.677, mud=0.525, sigma=120,
VEQ=1.22, Vstar=1e-6, Dc=1e-5, a=0.01, fa=0.1, b=0.015,
dtauf = -2.6,
ynuc=8000, znuc=-7500, wnuc=1500, segnuc=1) {
taup = mus * sigma[1]
tauf = mud * sigma[1]
aprime = fa*a
mu0 = (tauf + dtauf)/sigma[1] + (b-a)*log(VEQ/Vstar)
theta0 = (Dc/Vstar)^((b-a)/(b-a+aprime)) * ((l*b*sigma[1])/(lameMu*VEQ))^((aprime)/(b-a+aprime)) *
exp( (taup/sigma[1] - mu0)/(b - a + aprime) )
theta.nuc = RStheta(Vstar, mu0, taup, sigma[1], a, b, VEQ, Dc)
theta = rep(theta0, flt$np)
nuc = which(flt$iseg == segnuc & abs(flt$z - znuc) <= wnuc & abs(flt$y - ynuc) <= wnuc)
theta[nuc] = theta.nuc
write(theta, file=file.path(rundir, "initTheta.txt"))
return(mu0)
}
buildRSQSimInputFile =
function(rundir, run, a, b, Dc, mu0, ddotStar, VEQ, lameLambda, lameMu) {
unlink("eqs..out")
system2("/Users/dinger/RSQSim.github.master/runRSQSim", stdout=NULL, stderr=NULL)
eqs = readEqs("eqs..out", saveRData=FALSE)
unlink(c("eqs..out", "eqs..out.RData", ".Screen.Time.out", ".times_saved.out"))
params = eqs$params
rm(eqs)
params$A_1 = a
params$B_1 = b
params$Dc_1 = Dc
params$mu0_1 = mu0
params$ddotStar_1 = Vstar
params$alpha_1 = 0.25
params$ddotEQ_1 = VEQ
params$lameLambda = lameLambda
params$lameMu = lameMu
params$nEq = 1
params$faultFname = "faults.flt"
params$outFnameInfix = paste0("TPV15.run", run)
params$writeTau = 2
params$writeSigma = 2
params$initTauFname = "initTau.txt"
params$initSigmaFname = "initSigma.txt"
params$initThetaFname = "initTheta.txt"
writeRSQSimInputFile(params, file.path(rundir, paste0("TPV15.run", run, ".in")))
}
Note that with 1 km elements and set up as below, rupture would not propagate away from nucleation region. Currently using 400 m elements.
run = 1
rundir = file.path("TPV15", paste0("Run", run))
system(paste("mkdir -p", rundir))
l = 400
setoff = 0
ynuc = 8000
znuc = -7500
wnuc = 1500
segnuc = 1
taunuc = 81.6
mus = 0.677
mud = 0.525
initSigma = c(120, 120)
initTau = c(70, 78)
VEQ = 1.22 # note: usual default is 1.0, shear impedance in this case says 1.57
Vstar = 1e-6
Dc = 1e-5
a = 0.01
fa = 0.1
b = 0.015
dtauf = -2.6
eqs = NULL
flt = buildFault(l=l, setoff=setoff, rundir=rundir)
initStresses(flt, initSigma, initTau,
ynuc=ynuc, znuc=znuc, wnuc=wnuc, segnuc=segnuc, taunuc=taunuc)
mu0 = fricTheta0(flt, rundir, mus=mus, mud=mud, sigma=initSigma,
VEQ=VEQ, Vstar=Vstar, Dc=Dc, a=a, fa=fa, b=b,
dtauf = dtauf,
ynuc=ynuc, znuc=znuc, wnuc=wnuc, segnuc=segnuc)
buildRSQSimInputFile(rundir=rundir, run=run, a=a, b=b, Dc=Dc, mu0=mu0,
ddotStar=ddotStar, VEQ=VEQ,
lameLambda=lameLambda, lameMu=lameMu)
## [1] "Warning: /Users/dinger/Proposals/2019SCEC-branch/SCEC_CVS/faults.in doesn't exist"
here = getwd()
setwd(rundir)
unlink("eqs*out")
system2("mpirun",
args=c("-np", 4, "/Users/dinger/RSQSim.github.master/runRSQSim",
paste0("TPV15.run", run, ".in")),
stderr=paste0("run", run, ".stderr"),
stdout=paste0("run", run, ".stdout"))
eqs[[run]] = readEqs(paste0("eqs.TPV15.run", run, ".out"))
setwd(here)
hist(eqs[[1]]$taupList, breaks="FD")
c(mean(eqs[[1]]$taupList), mus*initSigma[1])
## [1] 81.47452 81.24000
tauf.rsq = readFaultVal(eqs[[1]], val="tau", patches=sort(eqs[[1]]$pList),
writeVal="final")$val[1,]
hist(tauf.rsq, breaks="FD")
c(mean(tauf.rsq), mud*initSigma[1])
## [1] 63.40268 63.00000
run = 2
rundir = file.path("TPV15", paste0("Run", run))
system(paste("mkdir -p", rundir))
setoff = 100
flt = buildFault(l=l, setoff=setoff, rundir=rundir)
initStresses(flt, initSigma, initTau,
ynuc=ynuc, znuc=znuc, wnuc=wnuc, segnuc=segnuc, taunuc=taunuc)
mu0 = fricTheta0(flt, rundir, mus=mus, mud=mud, sigma=initSigma,
VEQ=VEQ, Vstar=Vstar, Dc=Dc, a=a, fa=fa, b=b,
dtauf = dtauf,
ynuc=ynuc, znuc=znuc, wnuc=wnuc, segnuc=segnuc)
buildRSQSimInputFile(rundir=rundir, run=run, a=a, b=b, Dc=Dc, mu0=mu0,
ddotStar=ddotStar, VEQ=VEQ,
lameLambda=lameLambda, lameMu=lameMu)
## [1] "Warning: /Users/dinger/Proposals/2019SCEC-branch/SCEC_CVS/faults.in doesn't exist"
here = getwd()
setwd(rundir)
unlink("eqs*out")
system2("mpirun",
args=c("-np", 4, "/Users/dinger/RSQSim.github.master/runRSQSim",
paste0("TPV15.run", run, ".in")),
stderr=paste0("run", run, ".stderr"),
stdout=paste0("run", run, ".stdout"))
eqs[[run]] = readEqs(paste0("eqs.TPV15.run", run, ".out"))
setwd(here)
run = 3
rundir = file.path("TPV15", paste0("Run", run))
system(paste("mkdir -p", rundir))
setoff = 200
flt = buildFault(l=l, setoff=setoff, rundir=rundir)
initStresses(flt, initSigma, initTau,
ynuc=ynuc, znuc=znuc, wnuc=wnuc, segnuc=segnuc, taunuc=taunuc)
mu0 = fricTheta0(flt, rundir, mus=mus, mud=mud, sigma=initSigma,
VEQ=VEQ, Vstar=Vstar, Dc=Dc, a=a, fa=fa, b=b,
dtauf = dtauf,
ynuc=ynuc, znuc=znuc, wnuc=wnuc, segnuc=segnuc)
buildRSQSimInputFile(rundir=rundir, run=run, a=a, b=b, Dc=Dc, mu0=mu0,
ddotStar=ddotStar, VEQ=VEQ,
lameLambda=lameLambda, lameMu=lameMu)
## [1] "Warning: /Users/dinger/Proposals/2019SCEC-branch/SCEC_CVS/faults.in doesn't exist"
here = getwd()
setwd(rundir)
unlink("eqs*out")
system2("mpirun",
args=c("-np", 4, "/Users/dinger/RSQSim.github.master/runRSQSim",
paste0("TPV15.run", run, ".in")),
stderr=paste0("run", run, ".stderr"),
stdout=paste0("run", run, ".stdout"))
eqs[[run]] = readEqs(paste0("eqs.TPV15.run", run, ".out"))
setwd(here)
run = 4
rundir = file.path("TPV15", paste0("Run", run))
system(paste("mkdir -p", rundir))
setoff = 300
flt = buildFault(l=l, setoff=setoff, rundir=rundir)
initStresses(flt, initSigma, initTau,
ynuc=ynuc, znuc=znuc, wnuc=wnuc, segnuc=segnuc, taunuc=taunuc)
mu0 = fricTheta0(flt, rundir, mus=mus, mud=mud, sigma=initSigma,
VEQ=VEQ, Vstar=Vstar, Dc=Dc, a=a, fa=fa, b=b,
dtauf = dtauf,
ynuc=ynuc, znuc=znuc, wnuc=wnuc, segnuc=segnuc)
buildRSQSimInputFile(rundir=rundir, run=run, a=a, b=b, Dc=Dc, mu0=mu0,
ddotStar=ddotStar, VEQ=VEQ,
lameLambda=lameLambda, lameMu=lameMu)
## [1] "Warning: /Users/dinger/Proposals/2019SCEC-branch/SCEC_CVS/faults.in doesn't exist"
here = getwd()
setwd(rundir)
unlink("eqs*out")
system2("mpirun",
args=c("-np", 4, "/Users/dinger/RSQSim.github.master/runRSQSim",
paste0("TPV15.run", run, ".in")),
stderr=paste0("run", run, ".stderr"),
stdout=paste0("run", run, ".stdout"))
eqs[[run]] = readEqs(paste0("eqs.TPV15.run", run, ".out"))
setwd(here)
Only continued to this offset because this is one of our element sizes, which what TPV15 did - but with smaller elements
run = 5
rundir = file.path("TPV15", paste0("Run", run))
system(paste("mkdir -p", rundir))
setoff = 400
flt = buildFault(l=l, setoff=setoff, rundir=rundir)
initStresses(flt, initSigma, initTau,
ynuc=ynuc, znuc=znuc, wnuc=wnuc, segnuc=segnuc, taunuc=taunuc)
mu0 = fricTheta0(flt, rundir, mus=mus, mud=mud, sigma=initSigma,
VEQ=VEQ, Vstar=Vstar, Dc=Dc, a=a, fa=fa, b=b,
dtauf = dtauf,
ynuc=ynuc, znuc=znuc, wnuc=wnuc, segnuc=segnuc)
buildRSQSimInputFile(rundir=rundir, run=run, a=a, b=b, Dc=Dc, mu0=mu0,
ddotStar=ddotStar, VEQ=VEQ,
lameLambda=lameLambda, lameMu=lameMu)
## [1] "Warning: /Users/dinger/Proposals/2019SCEC-branch/SCEC_CVS/faults.in doesn't exist"
here = getwd()
setwd(rundir)
unlink("eqs*out")
system2("mpirun",
args=c("-np", 4, "/Users/dinger/RSQSim.github.master/runRSQSim",
paste0("TPV15.run", run, ".in")),
stderr=paste0("run", run, ".stderr"),
stdout=paste0("run", run, ".stdout"))
eqs[[run]] = readEqs(paste0("eqs.TPV15.run", run, ".out"))
setwd(here)
run = 6
rundir = file.path("TPV15", paste0("Run", run))
system(paste("mkdir -p", rundir))
l = 200
setoff = 0
VEQ = 1.22 # note: usual default is 1.0, shear impedance in this case says 1.57
b = 0.015
dtauf = -2.6
flt = buildFault(l=l, setoff=setoff, rundir=rundir)
initStresses(flt, initSigma, initTau,
ynuc=ynuc, znuc=znuc, wnuc=wnuc, segnuc=segnuc, taunuc=taunuc)
mu0 = fricTheta0(flt, rundir, mus=mus, mud=mud, sigma=initSigma,
VEQ=VEQ, Vstar=Vstar, Dc=Dc, a=a, fa=fa, b=b,
dtauf = dtauf,
ynuc=ynuc, znuc=znuc, wnuc=wnuc, segnuc=segnuc)
buildRSQSimInputFile(rundir=rundir, run=run, a=a, b=b, Dc=Dc, mu0=mu0,
ddotStar=ddotStar, VEQ=VEQ,
lameLambda=lameLambda, lameMu=lameMu)
## [1] "Warning: /Users/dinger/Proposals/2019SCEC-branch/SCEC_CVS/faults.in doesn't exist"
here = getwd()
setwd(rundir)
unlink("eqs*out")
system2("mpirun",
args=c("-np", 4, "/Users/dinger/RSQSim.github.master/runRSQSim",
paste0("TPV15.run", run, ".in")),
stderr=paste0("run", run, ".stderr"),
stdout=paste0("run", run, ".stdout"))
eqs[[run]] = readEqs(paste0("eqs.TPV15.run", run, ".out"))
setwd(here)
hist(eqs[[run]]$taupList, breaks="FD")
c(mean(eqs[[run]]$taupList), mus*initSigma[1])
## [1] 81.56892 81.24000
tauf.rsq = readFaultVal(eqs[[run]], val="tau", patches=sort(eqs[[run]]$pList),
writeVal="final")$val[1,]
hist(tauf.rsq, breaks="FD")
c(mean(tauf.rsq), mud*initSigma[1])
## [1] 63.7541 63.0000
run = 7
rundir = file.path("TPV15", paste0("Run", run))
system(paste("mkdir -p", rundir))
setoff = 100
flt = buildFault(l=l, setoff=setoff, rundir=rundir)
initStresses(flt, initSigma, initTau,
ynuc=ynuc, znuc=znuc, wnuc=wnuc, segnuc=segnuc, taunuc=taunuc)
mu0 = fricTheta0(flt, rundir, mus=mus, mud=mud, sigma=initSigma,
VEQ=VEQ, Vstar=Vstar, Dc=Dc, a=a, fa=fa, b=b,
dtauf = dtauf,
ynuc=ynuc, znuc=znuc, wnuc=wnuc, segnuc=segnuc)
buildRSQSimInputFile(rundir=rundir, run=run, a=a, b=b, Dc=Dc, mu0=mu0,
ddotStar=ddotStar, VEQ=VEQ,
lameLambda=lameLambda, lameMu=lameMu)
## [1] "Warning: /Users/dinger/Proposals/2019SCEC-branch/SCEC_CVS/faults.in doesn't exist"
here = getwd()
setwd(rundir)
unlink("eqs*out")
system2("mpirun",
args=c("-np", 4, "/Users/dinger/RSQSim.github.master/runRSQSim",
paste0("TPV15.run", run, ".in")),
stderr=paste0("run", run, ".stderr"),
stdout=paste0("run", run, ".stdout"))
eqs[[run]] = readEqs(paste0("eqs.TPV15.run", run, ".out"))
setwd(here)
run = 8
rundir = file.path("TPV15", paste0("Run", run))
system(paste("mkdir -p", rundir))
setoff = 200
flt = buildFault(l=l, setoff=setoff, rundir=rundir)
initStresses(flt, initSigma, initTau,
ynuc=ynuc, znuc=znuc, wnuc=wnuc, segnuc=segnuc, taunuc=taunuc)
mu0 = fricTheta0(flt, rundir, mus=mus, mud=mud, sigma=initSigma,
VEQ=VEQ, Vstar=Vstar, Dc=Dc, a=a, fa=fa, b=b,
dtauf = dtauf,
ynuc=ynuc, znuc=znuc, wnuc=wnuc, segnuc=segnuc)
buildRSQSimInputFile(rundir=rundir, run=run, a=a, b=b, Dc=Dc, mu0=mu0,
ddotStar=ddotStar, VEQ=VEQ,
lameLambda=lameLambda, lameMu=lameMu)
## [1] "Warning: /Users/dinger/Proposals/2019SCEC-branch/SCEC_CVS/faults.in doesn't exist"
here = getwd()
setwd(rundir)
unlink("eqs*out")
system2("mpirun",
args=c("-np", 4, "/Users/dinger/RSQSim.github.master/runRSQSim",
paste0("TPV15.run", run, ".in")),
stderr=paste0("run", run, ".stderr"),
stdout=paste0("run", run, ".stdout"))
eqs[[run]] = readEqs(paste0("eqs.TPV15.run", run, ".out"))
setwd(here)
run = 9
rundir = file.path("TPV15", paste0("Run", run))
system(paste("mkdir -p", rundir))
setoff = 300
flt = buildFault(l=l, setoff=setoff, rundir=rundir)
initStresses(flt, initSigma, initTau,
ynuc=ynuc, znuc=znuc, wnuc=wnuc, segnuc=segnuc, taunuc=taunuc)
mu0 = fricTheta0(flt, rundir, mus=mus, mud=mud, sigma=initSigma,
VEQ=VEQ, Vstar=Vstar, Dc=Dc, a=a, fa=fa, b=b,
dtauf = dtauf,
ynuc=ynuc, znuc=znuc, wnuc=wnuc, segnuc=segnuc)
buildRSQSimInputFile(rundir=rundir, run=run, a=a, b=b, Dc=Dc, mu0=mu0,
ddotStar=ddotStar, VEQ=VEQ,
lameLambda=lameLambda, lameMu=lameMu)
## [1] "Warning: /Users/dinger/Proposals/2019SCEC-branch/SCEC_CVS/faults.in doesn't exist"
here = getwd()
setwd(rundir)
unlink("eqs*out")
system2("mpirun",
args=c("-np", 4, "/Users/dinger/RSQSim.github.master/runRSQSim",
paste0("TPV15.run", run, ".in")),
stderr=paste0("run", run, ".stderr"),
stdout=paste0("run", run, ".stdout"))
eqs[[run]] = readEqs(paste0("eqs.TPV15.run", run, ".out"))
setwd(here)
run = 10
rundir = file.path("TPV15", paste0("Run", run))
system(paste("mkdir -p", rundir))
setoff = 400
flt = buildFault(l=l, setoff=setoff, rundir=rundir)
initStresses(flt, initSigma, initTau,
ynuc=ynuc, znuc=znuc, wnuc=wnuc, segnuc=segnuc, taunuc=taunuc)
mu0 = fricTheta0(flt, rundir, mus=mus, mud=mud, sigma=initSigma,
VEQ=VEQ, Vstar=Vstar, Dc=Dc, a=a, fa=fa, b=b,
dtauf = dtauf,
ynuc=ynuc, znuc=znuc, wnuc=wnuc, segnuc=segnuc)
buildRSQSimInputFile(rundir=rundir, run=run, a=a, b=b, Dc=Dc, mu0=mu0,
ddotStar=ddotStar, VEQ=VEQ,
lameLambda=lameLambda, lameMu=lameMu)
## [1] "Warning: /Users/dinger/Proposals/2019SCEC-branch/SCEC_CVS/faults.in doesn't exist"
here = getwd()
setwd(rundir)
unlink("eqs*out")
system2("mpirun",
args=c("-np", 4, "/Users/dinger/RSQSim.github.master/runRSQSim",
paste0("TPV15.run", run, ".in")),
stderr=paste0("run", run, ".stderr"),
stdout=paste0("run", run, ".stdout"))
eqs[[run]] = readEqs(paste0("eqs.TPV15.run", run, ".out"))
setwd(here)